#include <cstdlib>
#include <ctime>
#include <cmath>
#include <new>
#include <cstring>
#include <limits>
#include "header.h"

class RandomNumGenerator
{
protected:
	unsigned long x1n,x2n,x1n_1,x2n_1,x1n_2,x2n_2;
public:
	RandomNumGenerator() {
		setSeed();
	};
	~RandomNumGenerator() {
		;
	};
	inline void setSeed(unsigned long a1, unsigned long a2, unsigned long a3, unsigned long a4, unsigned long a5, unsigned long a6) {
		x1n_1=a1;
		x1n_2=a2;
		x1n=a3;
		x2n_1=a4;
		x2n_2=a5;
		x2n=a6;
	};
	inline void setSeed() {
		srand(time(NULL));
		x1n_1=rand();
		x1n_2=rand();
		x1n=rand();
		x2n_1=rand();
		x2n_2=rand();
		x2n=rand();
	};
	inline void setSeed(unsigned seed) {
		srand(seed);
		x1n_1=seed;
		x1n_2=rand();
		x1n=rand();
		x2n_1=rand();
		x2n_2=rand();
		x2n=rand();
	};
	inline double next(){
		unsigned long z,x1n_3,x2n_3;
		x1n_2=x1n_1;
		x2n_2=x2n_1;
		x1n_3=x1n_2;
		x2n_3=x2n_2;
		x1n_1=x1n;
		x2n_1=x2n;
		return ((z=((x1n=(1403580u*x1n_2-810728u*x1n_3)%4292967087u)-(x2n=(527612u*x2n_1-1370589u*x2n_3)%4294944443u))%4294967087u)==0)?(4294967087u/(double)4294967088u):(z/(double)4294967087u);
	};
	inline double unif(double a,double b) {
		return next()*(b-a)+a;
	};
};

RandomNumGenerator randNG;

double expoRand(double lambda) {
	return -log(1-randNG.next())/lambda;
}

double unifRand(double a, double b) {
	return randNG.unif(a,b);
}

void setSeed(unsigned seed) {
	randNG.setSeed(seed);
}

void PoissonProcess::doubleSize() {
	double *ptr;
	unsigned s;

	if ( mElem >= std::numeric_limits<unsigned>::max()/2 ) {
		throw LOGICERROR_TOOMANYEVENT;
	}

	try {
		ptr = new double[ mElem*2 ];
	}
	catch (std::bad_alloc &e) {
		throw MEMERROR_ALLOC;
	}
	if (ptr==NULL) {
		throw MEMERROR_ALLOC;
	}

	memcpy(ptr, eTime, mElem*sizeof(double));
	try {
		delete []eTime;
	}
	catch (...) {
		throw MEMERROR_FREE;
	}
	eTime=ptr;
	ptr=NULL;
	s=mElem;
	mElem*=2;
	while (s < mElem) {
		eTime[s]= eTime[s-1] + expoRand(lambda);
		s++;
	}
}
PoissonProcess::PoissonProcess (double lam) {
	unsigned char i;
	
	mElem=128;
	lambda=lam;
	
	try {
		eTime=(double *)new double[128];
	}
	catch (std::bad_alloc &e) {
		throw MEMERROR_ALLOC;
	}
	if (eTime==NULL) {
		throw MEMERROR_ALLOC;
	}

	eTime[0]=expoRand(lambda);

	for (i=1; i<mElem; i++) {
		eTime[i]=eTime[i-1]+expoRand(lambda);
	}
}
PoissonProcess::~PoissonProcess() {
	if (eTime!=NULL) {
		try {
			delete []eTime;
		}
		catch (...) {
			throw MEMERROR_FREE;
		}
	}
}
double PoissonProcess::getLambda() {
	return lambda;
};
unsigned PoissonProcess::value (double time) {
	if ( time < 0 ) {
		return 0u;
	}
	if ( time < eTime[mElem-1] ) {
		unsigned i,j;
		i=0;
		j=mElem-1;
		while ( j-i > 1 ) {
			unsigned m;
			m=(i+j)/2;
			if (eTime[m] > time) {
				j=m;
			}
			else {
				i=m;
			}
		}
		return i;
	}
	doubleSize();
	return value(time);
}
double PoissonProcess::event(unsigned n) {
	if ( n <= mElem ) {
		return eTime[n-1];
	}
	doubleSize();
	return event(n);
}
